This is an Electronic Supplement to the manuscript Marques et al. “Quantifying Deepwater Horizon oil spill induced injury on pelagic cetaceans” submitted to Marine Ecology Progress Series (MEPS).
There are 8 Electronic Supplements to the paper. The master file containing links to all the other 7 additional Electronic Supplements related to this paper is ES0_ElectronicSupplements.
You might be reading this file as a pdf or as an html. The links on this file only work if you are using the html version of it, available via the github repository or if you compiled it yourself as html and you have all the 8 html files in the same folder. Otherwise, as a pdf distributed as an Electronic Supplement to the MEPS paper, the links might not work. They might work. If it is possible, we can work with the MEPS Editorial Office such that we can add links below that will link to actual files, say the pdfs of each of these 8 files, on the publisher server.
This section details the version history for static pdf files submitted as Electronic Supplement pdfs:
1.0 [12 Aug 2022] Version included as a pdf Electronic Supplement in the MEPS original submission
2.0 [10 Feb 2023] Version included as a pdf Electronic Supplement in the MEPS re-submission after 1st round of reviewer’s comments
In this document we present the results of the age, sex and class structured model implemented for all species/stock in the GOM.
We report here the injury metrics originally considered in Schwacke et al. (2017): (1) lost cetacean years (LCY), the difference between the baseline and injured population sizes, summed over the entire modeled time period; (2) years to recovery (YTR), the number of years required before the injured population trajectory reaches 95% of the baseline population trajectory; and (3) maximum proportional decrease (MPD), the difference between the two population trajectories when the injured trajectory is at its lowest point, divided by the baseline. Note that LCY is intuitively the metric that is most dependent on initial population size.
The code below runs the simulations for all species, and stores the results in appropriate files, stored in appropriate species specific folders.
This code is included for completeness and does not need to be run as part of compiling this dynamic report as the results are already available in relevant files and it takes a considerably log time to run.
Creating objects to hold summary results:
The taxonomic units and the corresponding codes considered in this document are:
We re-order the results tables upfront so that results are organized by alphabetical order of the 4 letter code used to describe the taxonomic unit (this is also the order that the species are reported in tables 1 to 3 of the offshore paper).
The data for initial population sizes and proportion exposed are
imported from a txt file produced inside ES2_ElectronicSupplements.
We represent in turn, for each species, the simulated population trajectories under the oil spill and under the baseline scenarios, as well as the histograms of the distributions of the 3 injury metrics for each species. To each histogram four dashed vertical lines were added:
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 2549.139215 9.725500 6.375102
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 1429.418506 4.898900 5.185119
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 1901.525708 0.000000 1.801292
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 3071.523882 2.137900 4.168256
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 11369.564397 6.586900 5.658282
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 1057.863912 1.156400 3.538044
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 1249.503621 2.308200 4.194471
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 2449.136975 7.290900 5.735186
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 2434.417513 1.640900 4.006926
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 33809.632317 0.635700 3.474522
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 35190.477524 2.435900 4.191836
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 17714.767745 11.196700 8.899829
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 3715.207476 7.108100 5.883191
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 7509.270139 0.000000 1.296869
Population trajectories
Mean for each of the injury metrics
## LCY YTR MPD
## 1320.899319 1.592200 3.816632
In this section we present some summary analysis of the results across the 15 taxonomic units considered.
The population size and the proportion exposed are \(a priori\) expected to be key determinants of injury, since the population size will have a direct impact on LCY and the proportion exposed should have a direct impact on all 3 measures.
| Sp | Nstart | Pexp |
|---|---|---|
| Bwsp | 3098 | 0.140 |
| Fatt | 2152 | 0.145 |
| Ggri | 3063 | 0.197 |
| Gmac | 2065 | 0.120 |
| Kosp | 2322 | 0.197 |
| Pele | 5784 | 0.152 |
| Pmac | 2561 | 0.191 |
| Satt | 81233 | 0.155 |
| Sbre | 4867 | 0.148 |
| Scly | 9065 | 0.082 |
| Scoe | 5011 | 0.247 |
| Sfro | 48688 | 0.058 |
| Slon | 16501 | 0.377 |
| Ttro | 15791 | 0.206 |
| Ttrs | 64897 | 0.177 |
The proportion exposed against the initial population size is shown here for the 15 taxonomic units
In the following table we present the means of the 3 injury metrics:
| Sp | LCY | MPD | YTR |
|---|---|---|---|
| Bwsp | 1321 | -0.038 | 1.6 |
| Fatt | 1250 | -0.042 | 2.3 |
| Ggri | 2449 | -0.057 | 7.3 |
| Gmac | 1058 | -0.035 | 1.2 |
| Kosp | 1429 | -0.052 | 4.9 |
| Pele | 3072 | -0.042 | 2.1 |
| Pmac | 2549 | -0.064 | 9.7 |
| Satt | 33810 | -0.035 | 0.6 |
| Sbre | 2434 | -0.040 | 1.6 |
| Scly | 1902 | -0.018 | 0.0 |
| Scoe | 3715 | -0.059 | 7.1 |
| Sfro | 7509 | -0.013 | 0.0 |
| Slon | 17715 | -0.089 | 11.2 |
| Ttro | 11370 | -0.057 | 6.6 |
| Ttrs | 35190 | -0.042 | 2.4 |
In the following table we present the medians and the 95% confidence intervals for our current estimates of injury. These correspond to the results in table 3 in the main paper.
| Sp | LCYmed | LCYl | LCYu | MPDmed | MPDl | MPDu | YTRmed | YTRl | YTRu |
|---|---|---|---|---|---|---|---|---|---|
| Bwsp | 1207 | 445 | 2899 | 3.683 | 1.386 | 6.951 | 0 | 0 | 11 |
| Fatt | 1122 | 437 | 2809 | 4.186 | 1.939 | 6.642 | 0 | 0 | 13 |
| Ggri | 2264 | 978 | 5049 | 5.623 | 2.569 | 9.494 | 8 | 0 | 18 |
| Gmac | 964 | 385 | 2291 | 3.435 | 1.486 | 6.216 | 0 | 0 | 12 |
| Kosp | 1294 | 519 | 3112 | 5.114 | 2.315 | 8.492 | 5 | 0 | 14 |
| Pele | 2761 | 1081 | 6767 | 4.140 | 1.925 | 6.601 | 0 | 0 | 12 |
| Pmac | 2356 | 996 | 5240 | 6.344 | 2.672 | 10.189 | 11 | 0 | 21 |
| Satt | 31372 | 12884 | 67606 | 3.433 | 1.563 | 5.615 | 0 | 0 | 9 |
| Sbre | 2232 | 901 | 5207 | 3.995 | 1.836 | 6.349 | 0 | 0 | 11 |
| Scly | 1726 | 664 | 4067 | 1.734 | 0.747 | 3.261 | 0 | 0 | 0 |
| Scoe | 3387 | 1369 | 7897 | 5.673 | 2.479 | 10.426 | 8 | 0 | 17 |
| Sfro | 6961 | 2722 | 15509 | 1.264 | 0.542 | 2.280 | 0 | 0 | 0 |
| Slon | 15255 | 5302 | 44903 | 8.363 | 3.277 | 17.716 | 12 | 0 | 21 |
| Ttro | 10537 | 4597 | 23220 | 5.610 | 2.546 | 9.125 | 8 | 0 | 16 |
| Ttrs | 32584 | 13377 | 71967 | 4.152 | 1.859 | 6.758 | 0 | 0 | 13 |
These are the values that are used in table 3 of the paper, and the code to do so is in the .Rmd.
Representing the above information visually, focusing on each of the 3 injury metrics:
In this section we explore the results obtained in terms of the different injury metrics and how these injury metrics might be explained by the different input parameters. To do so we consider linear models to model the mean value of the injury per species as a function of the mean values of the input parameters. After a preliminary exploratory analysis, only the following four parameters seem to have a non-negligible potential relevant effect on the injury metrics:
Here we plot the values of the injury metrics as a function of the four key parameters going into the model.
The next three plots describe the relationship between the proportion exposed and Lost Cetacean Years, the Maximum Proportion Decrease and the Years to recovery, respectively
The next three plots describe the relationship between the initial population size and Lost Cetacean Years, the Maximum Proportion Decrease and the Years to recovery, respectively.
The next three plots describe the relationship between the gestation duration and Lost Cetacean Years, the Maximum Proportion Decrease and the Years to recovery, respectively
The next three plots describe the relationship between the survival reduction and Lost Cetacean Years, the Maximum Proportion Decrease and the Years to recovery, respectively
Here we consider standard linear models to explain the injury metrics as a function of the input parameters. To obtain what we considered the most parsimonious model we conducted a traditional model selection approach using AIC, selecting the model with the lowest AIC.
We considered only the 4 explanatory variables which, given the exploratory data analysis above, could explain the injury metrics. Given only 4 explanatory variables, we conducted a full search over the possible models, ignoring interactions, using function bestglm from the package bestglm (McLeod & Lai, 2020).
The summary of the best linear model follows.
## AIC
## BICq equivalent for q in (0.00257603089392322, 0.715027970358707)
## Best Model:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8585.8086023 2457.09226124 -3.494296 0.00442815340683
## N 0.4263961 0.03567448 11.952412 0.00000005054783
## p 54656.6056191 12278.22336920 4.451508 0.00079078410382
##
## Call:
## lm(formula = y ~ p + N, data = LCY4bestglm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7835.4 -1189.7 855.8 1477.7 6430.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8585.80860 2457.09226 -3.494 0.004428 **
## p 54656.60562 12278.22337 4.452 0.000791 ***
## N 0.42640 0.03567 11.952 0.0000000505 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3380 on 12 degrees of freedom
## Multiple R-squared: 0.9263, Adjusted R-squared: 0.9141
## F-statistic: 75.45 on 2 and 12 DF, p-value: 0.0000001598
We can look at the ANOVA table for the most parsimonious model based on lowest AIC criteria.
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## p 1 91857453 91857453 8.0395 0.01502 *
## N 1 1632291730 1632291730 142.8602 0.00000005055 ***
## Residuals 12 137109617 11425801
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Using the most parsimonious model, the proportion of the variability explained by each input retained in the regression model and the remaining unexplained variation can be obtained from the ANOVA table of the linear regression (as shown above). The proportion of variance explained in Lost Cetacean Years by the initial population size is 87.7 %, while the proportion exposed accounts for 4.94 %. Only 7.37 of the variability remains unexplained.
We illustrate how well the model predicts the observed values in the plot below:
The summary of the best linear model follows.
## AIC
## BICq equivalent for q in (0.0565670801885912, 0.714907105110615)
## Best Model:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -14.46579200 3.864204450 -3.743537 0.002803820209
## gd 0.02778531 0.009303748 2.986464 0.011348598337
## p 44.99620475 5.546928494 8.111914 0.000003260501
##
## Call:
## lm(formula = y ~ p + gd, data = YTR4bestglm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.364 -1.229 0.256 1.220 2.177
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -14.465792 3.864204 -3.744 0.0028 **
## p 44.996205 5.546928 8.112 0.00000326 ***
## gd 0.027785 0.009304 2.986 0.0113 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.515 on 12 degrees of freedom
## Multiple R-squared: 0.8502, Adjusted R-squared: 0.8253
## F-statistic: 34.06 on 2 and 12 DF, p-value: 0.00001129
We can look at the ANOVA table for the most parsimonious model based on lowest AIC criteria.
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## p 1 135.955 135.955 59.199 0.000005592 ***
## gd 1 20.483 20.483 8.919 0.01135 *
## Residuals 12 27.559 2.297
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The proportion of variance explained in Years to Recovery by the proportion exposed is 73.89 %, while the gestation duration accounts for 11.13 %. About 14.98 % of the total variability remains unexplained.
We illustrate how well the model predicts the observed values in the plot below:
The summary of the best linear model follows.
## AIC
## BICq equivalent for q in (0.000211081585492212, 0.79448839248248)
## Best Model:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.14127227832 0.03051310684 -4.629888 0.0007285396508132729
## gd -0.00006689489 0.00001223598 -5.467064 0.0001957224385135165
## sr 0.19516442031 0.03283771043 5.943302 0.0000968549623446209
## p -0.25059730715 0.00616856884 -40.624870 0.0000000000002441672
##
## Call:
## lm(formula = y ~ p + gd + sr, data = MPD4bestglm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0020501 -0.0011974 0.0001799 0.0007734 0.0033525
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.14127228 0.03051311 -4.630 0.000729 ***
## p -0.25059731 0.00616857 -40.625 0.000000000000244 ***
## gd -0.00006689 0.00001224 -5.467 0.000196 ***
## sr 0.19516442 0.03283771 5.943 0.000096854962345 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.001685 on 11 degrees of freedom
## Multiple R-squared: 0.9936, Adjusted R-squared: 0.9918
## F-statistic: 567.1 on 3 and 11 DF, p-value: 0.000000000002469
We can look at the ANOVA table for the most parsimonious model based on lowest AIC criteria.
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## p 1 0.0044337 0.0044337 1561.410 0.0000000000003305 ***
## gd 1 0.0002966 0.0002966 104.451 0.0000005944763062 ***
## sr 1 0.0001003 0.0001003 35.323 0.0000968549623446 ***
## Residuals 11 0.0000312 0.0000028
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The proportion of variance explained in MPD by the proportion exposed is 91.19, while the gestation duration and the survival reduction only account for only 6.1 % and 2.06 %, respectively. About 0.64 % of the variability remains unexplained.
We illustrate how well the model predicts the observed values in the plot below:
DWH MMIQT 2015, Models and analyses for the quantification of injury to Gulf of Mexico cetaceans from the Deepwater Horizon Oil Spill, MM_TR.01_Schwacke_Quantification.of.lnjury.to.GOM.Cetaceans LINK.
McLeod A, Xu C, Lai Y (2020). bestglm: Best Subset GLM and Regression Utilities. R package version 0.37.3, https://CRAN.R-project.org/package=bestglm.
Schwacke, L.H., L. Thomas, R.S. Wells, W.E. McFee, A.A. Hahn, K.D. Mullin, E.S. Zolman, B.M. Quigley, T.K. Rowles and J.H. Schwacke. 2017. An age-, sex- and class-structured population model for estimating nearshore common bottlenose dolphin injury following the Deepwater Horizon oil spill. Endangered Species Research 33: 265-279. DOI: 10.3354/esr00777.